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Abstract 

The DELPHI detector at LEP has been used to measure multi-muon bun- 
dles originating from cosmic ray interactions with air. The cosmic events were 
recorded in "parasitic mode" between individual e+e~ interactions and the to- 
tal live time of this data taking is equivalent to L6 • 10^ seconds. The DELPHI 
apparatus is located about 100 metres underground and the 84 metres rock 
overburden imposes a cut-off of about 52 GeV/c on muon momenta. The data 
from the large volume Hadron Calorimeter allowed the muon multiplicity of 
54201 events to be reconstructed. The resulting muon multiplicity distribu- 
tion is compared with the prediction of the Monte Carlo simulation based on 
CORSIKA/QGSJETOl. The model fails to describe the abundance of high 
multiplicity events. The impact of QGSJET internal parameters on the results 
is also studied. 
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1 Introduction 

The DELPHI (DEtector with Lepton Photon and Hadron Identification) at CERN 
LEP (Large Electron Positron colUder) measured cosmic muons regularly in order to 
align and calibrate various subdetectors. A major upgrade of the DELPHI hadron 
calorimeter was completed in 1997. As a result the calorimeter granularity increased 
substantially and spectacular events like the one shown in Fig. 1 were registered. The 
trigger studies performed during 1998 have shown that DELPHI can register cosmic 
events during regular data taking. Whenever there was no triggered e^e^ interaction, 
the detector stayed active to record possible cosmic events. In this regime we were able 
to collect data throughout the years 1999 and 2000. 

The experimental hall of DELPHI was located 100 metres underground and the 
overburden imposed a cut-off of 52 GeV/c on the momenta of vertical muons. This, 
depending on the particular interaction model, corresponds to a lower limit of primary 
particle energies of about 10^"^ eV. The upper limit of primary energy, less than 10^^ eV, 
follows from the total measurement time of 1.6 ■ 10^ seconds. Although this hve time is 
small compared to standard cosmic ray experiments, the granularity of the detector and 
the momentum cut-off make the data interesting. The high energy muons originate from 
meson decays and other processes which take place in the upper atmosphere. They carry 
information about the first stages of the shower development. Consequently, these data 
reflect different aspects of the shower than those recorded by experiments on the ground, 
where the vast majority of detected muons originates from pion decays at low energies. 

Reconstruction of cosmic ray interactions at very high energies relics heavily on Monte 
Carlo (MC) simulations. Hence the interpretation of measured data is dependent on the 
models of shower propagation, including simulations of high energy hadron collisions, 
hadron decays and further development of the electromagnetic and hadronic compo- 
nents. While the particle decays and the shower propagation are well described, the 
most important source of uncertainties originates from models describing the high energy 
interactions of hadrons at the beginning of shower development. The interaction mod- 
els such as NEXUS [1], QGSJET [2] or SIBYLL [3] are tuned to available accelerator 
data at lower energies than those discussed in this paper. The collider experiments are 
more suited to study phenomena at larger transverse momenta. Thus our data, which 
can reveal features of particle interactions in the very forward region, are in this sense 
complementary. 

The muon component of cosmic ray showers has been studied with large ground ar- 
rays (e.g. [4,5]) or at large depths corresponding to a momentum cut-off above 1 TeV 
(e.g. [6-8]). The data at intermediate depths underground are scarce and the experi- 
ments detecting muons with a momentum cut-off around 100 GeV/c (e.g. [9]) use less 
precise detectors than the LEP experiments. Besides DELPHI, similar studies of cos- 
mic rays were performed at ALEPH [10] and L3+C [11]. Detailed model tests [12] show 
that QGSJET describes best the various correlations between hadronic, electromagnetic 
and muon components of atmospheric showers in the case of ground experiments. Data 
registered by underground experiments reflect different shower properties. The aim of 
this work is to test the interaction model, which sufficiently well describes the ground 
measurements, using multi-muon data detected underground. 

The detector and its overburden are described in Section 2. The conditions of event 
registration are mentioned in Section 3 and the procedure of event reconstruction is 
described in Section 4. The chain of programs used to simulate showers is described in 
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Section 5. The results obtained are given together with predictions of hadronic interaction 
models in Section 6 and they are discussed in the final Section 7. 

2 Detector and its location 

DELPHI was a classical collider experiment with numerous subdctcctors and a 
solenoidal magnetic field. A detailed description of the apparatus can be found in [13]. 
Only a few subdetectors were used for the cosmic muon detection, namely: Time Pro- 
jection Chamber (TPC), Time Of Fhght scintillation detector (TOP), Outer Detector 
(OD), Barrel part of HAdron calorimeter (HAB) and Barrel MUon chambers (MUB). All 
these parts were located in the barrel part of the detector (Fig. 2). TOP served to trigger 
cosmic events. 

The HAB detector was a sampling calorimeter and it contained 12000 limited streamer 
tubes. The iron of the magnet yoke served as an absorber. It consisted of 20 slabs 5 cm 
thick. Streamer tubes were inserted into the 2 cm wide gaps between individual iron 
plates. The gas mixture inside the tubes was composed of Ar(10%), CO2(60%) and iso- 
butane(30%). HAB with its large volume served as the backbone of muon detection. 
The detection area of HAB was 75 m^ in the horizontal plane. Each tube in the barrel 
part of the hadron calorimeter had an effective length of 3.6 m and its cross-section was 
1x8 cm^. All the tubes were parallel to the beam pipe. During the upgrade of the 
hadron calorimeter in the years 1995 - 1997 each tube was equipped with read-out of its 
cathode, which consisted of resistive varnish of the whole tube interior [14] . The smallest 
sensitive cell before the upgrade was about 20 x 30 x 35 cm^ in (9, 0, R) standard DELPHI 
coordinate system^ and the cells were organised in towers pointing to the centre of the 
detector. After the upgrade the cell size of the cathode readout in the barrel became 
360 X 8 X 7 cm^ [15, 16]. Consequently the granularity in the plane perpendicular to the 
beams increased about 14 times. Due to technical hmitations it was possible to read 
out signals only on the two outer front-ends of the barrel. The charge deposited on the 
cathode was integrated for 350 and accepted or rejected by a discriminator. Thus in 
this system of cathode read-out, the signals from individual tubes were either yes or no 
and the reconstructed tracks are in fact only projections of the muon trajectories onto 
the plane perpendicular to the LEP beams, separately for each half of HAB. 

The TPC was able to measure the full direction of muon tracks. Due to its relatively 
small volume it contained only a small fraction of the muons passing through DELPHI 
(TPC had 10 times smaller detection area compared to HAB). During the standard 
recording of e'^e~ collisions, the drift time in the TPC is measured from to which is given 
by the instant of beam cross-over (BCO) inside DELPHI. In the case of cosmic events 
to was the average arrival time of tracks to the OD. 

In extreme cases 50% or more of the tubes in one or both sides of HAB were hit. This 
led to saturated events where counting of individual muons was not possible anymore. 
However, the cosmic origin of these events is guaranteed, because in this case vacant 
tubes appear in parallel lines which follow the direction of the muon bundle and they 
cannot be caused by any noise in HAB. Moreover, in a few such events the lower bound 
on the number of muons could be roughly assessed from MUB. 

The apparatus was situated about 100 m underground. The surface altitude was 428 m 
above the sea level. The composition of the rock above the DELPHI experiment is known 
from a geological survey performed for civil engineering purposes. The simplified picture 

^as defined e.g. in [13] - R radius, (p azimuth angle in plane perpendicular to the beam pipe and 6 polar angle (= 
along beam) 
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of the overburden structure could be approximated by 5 major geological layers with 
different mass densities. The density of the rock in the vertical direction varies between 
2.2 g/cm^ and 2.5 g/cm^ depending on the layer. The total vertical depth of DELPHI 
location is about 19640 g/cvo?. The resulting energy cut-off for vertical cosmic muons is 
~ 52 GeV. The detector was located in a large experimental cavern equipped with three 
access shafts shown in Fig. 3. This scheme of the experimental area and the overburden 
was used in simulations. 

3 Trigger 

The trigger of cosmic events was entirely based on TOF. This detector consisted of a 
single layer of plastic scintillation counters. Each one was read out by two photomultipli- 
ers. The scintillator planks covered the internal side of HAB. Initial attempts to trigger 
on single muons led to a high trigger rate. Therefore in 1999 the trigger was set up to 
demand at least 3 active detector sectors to accept an event. It ran in so-called "parasitic 
mode", i.e. whenever there was no triggered e^e~ interaction, the trigger stayed sensitive 
to cosmic events for 4.1 fis after each beam crossing. This short detection window was 
optimised for e^e~ interactions. 

The beam crossing frequency depended on the number of e^{e^) bunches in the col- 
lider. During the running mode with 4 bunches in the machine, the beam crossing period 
was ~ 22.2 yus, while in the 8 bunch mode the period decreased to ~ 11.1 ^s. Con- 
sequently, the detector was sensitive to cosmic events for 18% of the total data taking 
time in 4 bunch mode and for 37% in 8 bunch mode. Dedicated cosmic runs (without 
the beams in the collider) have been performed mainly at the beginning of each year. 
Although there were no e+e" collisions, BCO signals were issued to mimic the 8 bunch 
mode. 

In an ideal case, two muons passing TOF would be sufficient to activate the trigger. 
In reality the TOF detection efficiency in Z° — > jJif^n' events was 84%. However, with 
increasing muon multiplicity the TOF trigger efficiency quickly approaches almost 100%. 
Already for muon multiplicity A^"^ = 5 the TOF efficiency is 99%, for lower multiplicities 
= 3(4) the corresponding efficiencies arc 94(97)%. It was found in [17] that with 5 
or more muons the trigger stability is assured. Fig. 4 plots the rate of events with muon 
multiplicity higher than 5 in different run periods. The event rates are consistent within 
statistical errors and there is no difference between the runs with and without beams 
in LEP. In total, taking into account various bunch schemes and the 4.1 fis detection 
window, the accumulated effective live time is Tgff = 1.6 ■ 10^ s (= 18.5 days). 

4 Event reconstruction 

The tracks of cosmic muons were reconstructed from hadron calorimeter data by the 
ECTANA program [18], which scans signals in the HAB modules and finds track patterns 
of hit streamer tubes. This package has the advantage that it was developed not only for 
studies of e'^e~ collisions, i.e. tracks coming from the interaction point in the centre of the 
detector, but it has the option for cosmic events as well. When running in cosmic mode 
it allows tracks originating anywhere in the calorimeter to be reconstructed without an 
explicit cut on the track impact parameter. The search for active streamer tubes starts 
from the outer planes of a given module and continues inwards. A group of at least 4 
aligned hits is taken as a track element. The track element is also required to have a 
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reasonable density of hits, at least 30% of tubes along its length have to be active. All 
possible hypotheses starting from a certain hit found during the scan are analysed, and 
the positions of hits are fitted by a straight hne. The best fit in terms of the number of 
hits and is stored. Before accepting the track, its similarity with other hypotheses was 
checked to avoid double counting. 

The length of the reconstructed track was required to be larger than 50 cm. It was 
possible to fit radii of curvature of the bent tracks, however, there were only a few such 
tracks and their radii were quite large. Therefore the coordinates of active tubes were fit- 
ted only by straight lines in the final analysis. The matching between track elements from 
different calorimeter sectors was performed. The number of reconstructed tracks was con- 
sidered as the reconstructed multiplicity of an event. The performance and functionality 
of the ECTANA program were checked with MC studies that compared parameters of 
reconstructed and injected events. However, no MC tuning of the reconstruction software 
was needed. 

The analysed data sample consists of 54201 events with muon multiplicities bigger 
than 3. They were registered during the years 1999 and 2000. The number of events 
with multiplicity above a given value is given in Tab. 1 and the differential multiplicity 
distribution is shown in Fig. 5. 

Altogether there were only 7 saturated events like the one depicted in Fig. 6 where 
more then 50% of the tubes were hit. In the case of saturated events vacant tubes make 
parallel line patterns which cannot result from a ghtch of the electronics. The saturated 
events are expected to have multiplicity higher than the highest multiplicity reconstructed 
from unsaturated events (A*"^ > 127). Moreover, in two of these events we were able to 
assess the lower limit of the multiplicity from the proportionality between the number 
of MUB anode hits and reconstructed muon multiplicity from HAB (Fig. 7). However, 
this procedure was not possible in all events. The MUB time window is only 5.9 /xs after 
BCO and the events coming at the end of trigger time window 4.1 /xs after BCO are not 
registered properly as the necessary drift time is 2.5 /is. 

In general, the muon tracks inside bundles are almost parallel as demonstrated in 
Fig. 8. In this picture we plot the angle a between the vertical direction and the track 
projection onto the plane perpendicular to the LEP beams. The track collinearity helped 
to find high multiplicity events originating from muon interactions close to the detector. 
The manual scanning was done on all events with > 30. Altogether we have rejected 
14 events with diverse directions of tracks. They correspond to 1.3% of the 1065 scanned 
events. The parallelism of reconstructed tracks was checked also by the cut that requires 
more than 50% of reconstructed tracks to be aligned within 5° of the mean value of all 
track angles in the event. This cut rejected the same events as the scanning procedure. 

As already mentioned above, the cathode read-out could not detect how many muons 
hit one single tube. Therefore at higher multiplicities muons start to shadow each other 
and the reconstructed multiplicity is in fact a lower limit of the real event multiplicity. 
However, even the highest reconstructed multiplicities around 120 are still strongly cor- 
related with the initial multiplicity as can be seen from Fig. 9, where the reconstructed 
multiplicity in MC data is plotted as a function of the number of muons injected into 
HAB. 

Unlike the hadron calorimeter, the TPC gives full spatial information on traversing 
muons. The drawback is its relatively small size. The track reconstruction from the 
TPC was possible with standard DELPHI software tools with the provision for start 
of the drift time (see Section 2). Due to the disproportion of TPC and HAB sizes, the 
respective multiplicities do not correlate well. However, we were able to reconstruct the 
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muon bundle directions from the TPC and to compare the multiphcities from the TPC 
with MC predictions [17]. 

5 Simulation 

To simulate the response of DELPHI to cosmic-ray induced showers, we have set 
up a chain of simulation programs. The high energy interactions were modelled by the 
QGSJETOl [2] program implemented within the CORSIKA [19] package^ . The rock 
above the DELPHI detector and the shape of the experimental cavern as well as the basic 
structures such as concrete walls and the three access shafts were represented according 
to Fig. 3 and simulated by GEANT3 [20]. Full simulation of the detector response was 
provided by the DELSIM [21] simulation package. 

As the chemical composition of cosmic rays is not well known, we have used only 
two limiting cases of hadron primary particles - protons and iron nuclei. Data sets 
were generated for both types of primary particles in 12 energy intervals 10^^ - S-IO^^ 
3-10^^ - 10^'^ cV, etc. up to 3-10^'' - 10^^ eV. The lowest energy interval barely contributes 
due to the muon energy cut-off of 52 GeV and the condition > 3. Also the highest 
energy interval contributes very little, if at all, because of the relatively short observation 
time. The lower energy limit depends on the interaction model and on the thickness of 
the overburden while the upper limit is given by the flux value used for normalisation 
and the observation time. As these two limits are not given reliably we have used a wider 
energy range for the simulations. 

All CORSIKA simulations were done without "thinning". At high energies {E > 
10^^ eV) the thinning option speeds up simulations of showers with billions of secondary 
particles by discarding a defined fraction of the secondaries and by ascribing the remaining 
particles certain weights. However, this option might introduce additional systematic 
errors. For this reason full event simulation was used in the analysis. 

The data samples were generated according to an energy dependence ~ E~'^ using 
the spectral index 7 = 1 in order to obtain sufficient representation of events at the 
upper part of the energy spectrum. Events were then re-weighted according to one of the 
assumed energy spectra (see below). 

Shower centres were smeared uniformly over a circular area with radius i? = 200 m 
around the DELPHI detector. This radius value was chosen as optimal because smaller 
R values led to an increased fraction of lost events with small muon multiplicities while 
larger radii would imply the necessity of using large data samples to produce enough 
events with high muon multiplicities. This is demonstrated in Fig. 10 which shows the 
stability of the simulated multiplicity distribution as a function of R. For each radius the 
ratio of occupancies in two adjacent bins in the final integral multiplicity^ distribution is 
plotted. With increasing values of R, the simulated multiplicity distribution stabilises. 
At R — 200 m the stability is reached at all simulated energies. Furthermore, the radius 
of 200 m ensures that the fraction of lost events at the lowest multiplicity A^^ = 4 is 
smaller than 0.5%. 

During the smearing of showers with E < 10^^ eV each shower was used 10 times. 
For higher energies the number of moves is 100. Taking 100 moves at energy > 10^® eV, 
one CORSIKA generated shower contributed to the simulated spectrum at > 45 on 

*^ First analyses with QGSJET model were performed with CORSIKA ver. 6.014 from March 2002. Later studies of 
QGSJET with modified parameters used CORSIKA ver. 6.031 from February 2004. It was checked that the results of 
the two simulations were independent of the CORSIKA version. 

^Defined in Section 6. 
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average only once. Since the events with A'^ > 45 are dominated by primary energies 
higher than 10^^ eV, the relatively high number of moves is, in fact, chosen optimally. 

The generated data set at A^^ > 45 (which corresponds roughly to £■ ~ 10^^ eV) 
was about 20 times larger than the real data sample. At lower multiplicities (i.e. lower 
energies) the samples were about equal. The stability of the results was also checked by 
doubling the size of the MC data sets. 

The normalisation of the simulated multiplicity distributions depends necessarily on 
the assumed energy spectrum of primary particles. Four spectra corresponding to dif- 
ferent lines in Fig. 11 were assumed. Lines 1, 2 and 3 all represent power law indices 
7 = 2.7 below the knee (-Ekncc = 3 • 10^^ eV) and 7 = 3.0 above the knee, thus they have 
the same shape of energy dependence and they differ by the total flux only. Assumption 
lb is defined by exponents 7 = 2.6 below and 7 = 3.05 above the knee. These spectral 
indices were used for tests of QGSJETOl with changed internal parameters. 

The most notable contributions to the systematic errors arc our imperfect knowledge 
of the overburden and due to a hardware effect which in certain situations caused cross- 
talk of the cathode read-out and appeared as a wider muon track that can shadow more 
muons than the normal track. The effect of inaccurate knowledge of the overburden 
was taken into account by changing the rock density by ±5% in all geological layers. 
Changes of multiplicity distribution induced by this density variations stay within 5%. 
The cross-talk has been studied in detail in Z° Ai^A*~ interactions. Based on this it was 
incorporated into the MC. The systematic error induced by this effect was checked in MC 
by using two options: one with full cross-talk simulation taken into account and another 
with this simulation switched off. It was found that the impact of cross-talk on the final 
multiplicity distribution is less than 5% of the number of events at high multiplicities. 
The upper bound of the possible live time error was estimated using the knowledge of 
the DELPHI dead time and it is about 2%. Due to the DELPHI magnetic field, another 
effect which might induce systematic error is the possible track matching inefficiency in 
the upper and the lower part of HAB for low energy muons. Assuming only straight 
fines in track reconstruction we could double count curved tracks. The effect was studied 
using the option of the ECTANA package that enabled to search also for curved tracks. It 
was found that the maximal impact on the final multiplicity distribution decreases with 
increasing multiplicity and it is about 8% for multiplicities below 15, 4% at integrated 
multiplicities larger than 20, 3% for multiplicities larger than 45 and 2% for multiplicities 
larger than 70. 

The overall systematic error is < 8% at high multiplicities (A^^ > 45) which is below 
the statistical uncertainty. More detailed discussion of the whole simulation is provided 
in [17]. 

6 Results 

6.1 Directions of muon bundles 

The most straightforward and MC-independent results are those concerning the di- 
rections of muon bundles. As explained already above, it was possible to reconstruct 
the full spatial direction of the tracks only from TPC data. As the TPC reconstruction 
depends on the mean arrival time to OD, we have selected higher multiplicity events 
with more than 15 muons in HAB and at least 4 reconstructed tracks in TPC. This cut 
corresponds to primary energies of about ~ 10^^ eV. The sky plot of event directions in 
galactic coordinates is shown in Fig. 12. The event direction is given as a mean direction 
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of individual muons and the pointing precision is a few degrees due to multiple scattering 
in the overburden, detector precision and unknown core position of the shower. There is 
no apparent clustering of events. 

The absence of point sources is demonstrated also by the dependence of the event rate 
on sidereal time. Fig. 13 shows no significant modulation of the rate during the sidereal 
day. The small dip disappears at higher multiplicities. 

The lack of point like anisotropics in the data justifies the assumption of uniform 
distribution of cosmic ray directions which is used in MC simulations. 

6.2 Muon multiplicities 

The shadowing effect reduces the number of reconstructed tracks when compared to 
the number of muons entering the calorimeter. In fact wc measure only a lower limit 
of the event multiplicity and therefore we plot the integrated multiplicity distributions 
where all events with given multiplicity or higher contribute to the corresponding bin. 
The measured distribution is plotted in Fig. 14a together with MC simulations of proton 
and iron induced showers. 

Taking into account that the composition of primary cosmic rays is light at energies 
~ 10^^ eV the data should follow the MC prediction for proton primary particles at small 
multiphcities. This behaviour is guaranteed only by fiux value 1 from Fig. 11. However, 
this value represents the upper limit of measured fiuxes. Taking into account the spread 
of fiux 1 - 3, we obtain for the MC prediction the bands demonstrated in Fig. 14b. 

Evidently even the highest fiux value combined with the assumption of pure iron 
primaries is not sufficient to describe the surplus of high multiplicity events. The excess 
of events in the region > 80 is 1.9 a (based on statistical errors) for fiux 1; assuming 
a more realistic fiux value 2, the discrepancy reaches about 3 a. One is tempted to 
interpret Fig. 14 as a convolution of proton and iron induced showers. However, this 
would mean that the primary particles at lower energies would be only protons while 
at the higher energies the primaries would be entirely iron nuclei. The contributions of 
individual energy bins in the case of iron primaries are detailed in Fig. 15. Fig. 16 shows 
the distribution of projected angle a measured in HAB as compared to MC event samples 
with A^"^ > 4 and > 20 respectively. The lower multiplicity corresponds to the point 
in Fig. 14 where data can be described by proton primaries. The second multiplicity 
interval represents the region where MC simulation of iron nuclei best describes the data. 

The saturated events appear in the simulation in the same way as in the data as 
events with more than 50% of the tubes hit. In the case of primary protons and fiux 1 
the number of MC saturated events is 1.1 ± 0.4. In the case of iron primaries the total 
number of expected saturated events is 3.3 ± 1.1 compared to 7 saturated events in the 
real data. 

Although we have tested only the QGSJET model, it is clear that the use of other 
models would lead to an even greater discrepancy as QGSJET predicts higher muon den- 
sities close to the shower core than other models do (e.g. SIBYLL or DPMJET [24]). 
Because of this, it was suggested [25] to test the sensitivity of the produced multiplic- 
ity spectra to QGSJET internal parameters. In [26] a set of QGSJETOl parameters is 
modified; namely the inelastic cross-section of p-p (p-N) is reduced and the elasticity of 
the collisions is increased. It is argued that such modifications can improve consistency 
between measurements of cosmic ray composition by experiments based on shower ar- 
rays and by Cerenkov or fiuorescence telescopes. Reference [26] suggests several possible 
modifications. In the following we will keep its notation and denote the tested model as 
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modification 3a. The result obtained with the modified QGSJET is compared with the 
data and with the original QGSJETOl in Fig. 17a. 

The model 3a enlarges the region where the data are between the predictions for proton 
and iron induced showers. In the case of unmodified QGSJETOl, the data reach the iron 
curve at multiplicity ~ 20. Using 3a, the data are consistent with a mixture of light 
and heavy components up to a multiplicity ~ 70. The slight event excess in data is still 
apparent at the highest multiplicities, however, now with somewhat smaller significance. 

The number of events at low muon multiplicities in the case of proton primaries and 
model 3a ( Fig. 17a) is now larger than in the data. The smaller and more realistic flux 
lb predicts a number of low multiplicity events consistent with the data as seen from 
Fig. 17b. At high multiplicities the model 3a predicts of course less events with flux lb 
than with flux assumption 1. However, the prediction of model 3a with spectrum lb is 
still above the prediction of QGSJET with flux 1. 

7 Conclusions 

The fine granularity hadron calorimeter of the DELPHI experiment was used to mea- 
sure multi-muon events originating from cosmic ray showers. The multiplicity distribution 
of muon bundles cannot be described by current Monte Carlo models in a satisfactory 
way. It is difficult to express the disagreement quantitatively as we have to use flux 
values measured elsewhere and also the chemical composition of initial particles is not 
well known. However, even the combination of extreme assumptions of highest measured 
flux value and pure iron spectrum fails to describe the abundance of high multiplicity 
events. Similar qualitative conclusions can be drawn from measurements of ALEPH [10] 
and L3+C [27], where muon bundles (up to multiplicity of about ~ 30) were studied in 
coincidence with the ground array signals. 

The tested QGSJET-based model with modifled cross-sections [26] performs somewhat 
better but it uses a value of the p-p total cross-section at the lowest limit allowed by 
CDF [28], E710 [29] and E811 [30] measurements. Justification of this assumption can 
be given only by future experiments. Hadron interactions at energies beyond the reach 
of accelerators are not very well known. Recently, also a more exotic explanation [31], 
based on the assumption of the presence of strangelets in cosmic rays, has been suggested 
to describe enhanced production of high multiplicity multi-muon events. 

The main conclusion is that the multi-muon data from cosmic ray showers detected 
at intermediate depths are quite sensitive to the dynamics of initial high energy interac- 
tions. In our case the primary collisions leading to high multiplicity events {Nfj, > 45) 
correspond to interactions at energies equivalent to about 5 TeV in the pp centre-of-mass 
system. This energy region has been so far inaccessible to laboratory measurements. 
However, even after LHC data become available, muon underground measurements have 
the potential to reveal some details of interactions in the very forward direction which 
are inaccessible to collider experiments. Thus they are important for the tuning of high 
energy interaction models which are indispensable for measurements and energy recon- 
struction of cosmic rays at even higher energies of the order 10^° eV, inaccessible to 
present and near future accelerators. 
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number of events 


A^,.>3 


54201 


> 30 


1065 


> 70 


78 


iV^ > 100 


24 



Table 1: Multiphcities of reconstructed events. 
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Figure 1: High multiplicity cosmic event as seen by hadron calorimeter. The number of 
reconstructed tracks was 127. 




Figure 2: The layout of the DELPHI detector; the hatched area represents the hadron 
calorimeter. Subdetectors used in this work are marked by larger letters. 
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Access shafts: 
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Figure 3: Schematic picture of rock overburden above DELPHI detector. 
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Figure 4: Event rates {N^ > 5) for different run periods. 
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Figure 5: Differential muon multiplicity distribution. 





Figure 6: A saturated event in the hadron calorimeter. Vacant tubes show voids in 
muon bundle. 
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5ure 7: Multiplicity reconstruction of saturated events from MUB data. 




Figure 8: Projected angle distribution in a high multiplicity event. 
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Figure 9: The correlation between injected and reconstructed numbers of muons in MC 
simulation. 
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Figure 10: Ratio of two adjacent bins (see legends inside the plots) of integral multiplicity 
distribution as a function of the parameter R. The plots correspond to iron induced 
vertical showers at a primary energy of 10^^ eV (a) and 10^^ eV (b). 
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Figure 11: Assumed fluxes compared to various measurements. The picture is taken 
from [22] and modified. The squares close to line 1 correspond to results of Hav- 
erah Park taken from [23]. The data points were added using the macro available at 
http://astroparticle.uchicago.edu/announce.html. Fluxes are multiplied by E'^-'^ . 
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Figure 12: Galactic coordinates of events with more than 15 tracks in HAB and more 
than 3 reconstructed tracks in TPC 
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Figure 13: The event rate versus the sidereal time expressed in degrees. Events with 
more than 5 reconstructed muons are taken into account. 
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Figure 14: Integrated multiphcity measured in HAB together with the result of the MC 
simulation of iron and proton induced showers with assumed flux 1 (a) and fluxes 1-3 
(b). 
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Figure 15: Contributions of different energy intervals to tfie final integral multiplicity 
distribution. Primary particles are iron nuclei. 
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Figure 16: Cosine of the projected angle a at A*"^ > 4 (a) and iV^ > 20 (b) for iron 
simulation (squares), data (full line) and proton simulation (diamonds). Normalisation 
of MC curves is done according to flux 1 from Fig. 11. 
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Figure 17: (a) The integral multiplicity distribution for QGSJET and modification 3a 
compared to data. Flux 1 is assumed, (b) The integral multiplicity distribution for the 
modification 3a compared to data. Flux Ihis assumed. 



